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Every sixth death in industrialized countries occurs because of cardiac arrhythmias like ventricular 
tachycardia (VT) and ventricular fibrillation (VP) . There is growing consensus that VT is associated 
with an unbroken spiral wave of electrical activation on cardiac tissue but VF with broken waves, 
spiral turbulence, spatiotemporal chaos and rapid, irregular activation. Thus spiral-wave activity 
in cardiac tissue has been studied extensively. Nevertheless, many aspects of such spiral dynamics 
remain elusive because of the intrinsically high-dimensional nature of the cardiac-dynamical system. 
In particular, the role of tissue heterogeneities in the stability of cardiac spiral waves is still being 
investigated. Experiments with conduction inhomogeneities in cardiac tissue yield a variety of 
results: some suggest that conduction inhomogeneities can eliminate VF partially or completely, 
leading to VT or quiescence, but others show that VF is unaffected by obstacles. We propose 
theoretically that this variety of results is a natural manifestation of a complex, fractal-like boundary 
that must separate the basins of the attractors associated, respectively, with spiral breakup and 
single spiral wave. We substantiate this with extensive numerical studies of Panfilov and Luo-Rudy 
I models, where we show that the suppression of spiral breakup depends sensitively on the position, 
size, and nature of the inhomogeneity. 

PACS numbers: 



I. INTRODUCTION 



The challenge of understanding the dynamics of spi- 
ral waves in excitable media is especially important in 
cardiac tissue where such waves are implicated in life- 
threatening arrhythmias such as ventricular tachycardia 
(VT) and ventricular fibrillation (VF)fll, U, i, i, H H, 
Anatomical reentry because of conduction inhomo- 
geneities in cardiac tissue, and functional reentry^, 
which result from wave propagation around transiently 
inexcitable regions, are crucial for the initiation of RS (a 
single rotating spiral wave) and the initiation and mainte- 
nance of ST (spiral turbulence with broken waves). But 
the precise ways in which spiral waves are affected by 
obstacles in ventricular tissue is still not clear [sj. Spiral 
waves form when waves of excitation circulate around an 
anatomical obstacle _9j. However, AUesie et al have 
shown that spiral-wave formation can also occur with a 
functionally determined heterogeneity in the tissue. The 
interaction of such a wave with an anatomical obstacle 
can be quite complex especially in the spatiotemporally 
chaotic state associated with spiral turbulence. Indeed, 
experiments with obstacles in cardiac tissue have yielded 
a variety of results. For example, some experiments 11 
report that small obstacles do not affect spiral waves but, 
as the size of the obstacle is increased, such a wave can 
get pinned to the obstacle. Various other experiments 
have discussed the role of an anatomical obstacle as an 
anchoring site for spiral waves, which can lead to the 
conversion of ST into RS [12, ^1^, [l3|. Davidcnko et al 
have found that, when they induced spiral waves 
in cardiac tissue preparations "... in most episodes, the 



spiral was anchored to small arteries or bands of con- 
nective tissue, and gave rise to stationary rotations. In 
some cases the core drifted away from its site of origin 
and dissipated at the tissue border." Other studies have 
shown [la EM 13 that an obstacle, in the path of a 
moving spiral wave, can break it and lead to many com- 
peting spiral waves. Recent experiments by Hwang et 
al [20| have suggested that multistability of spirals with 
different periods in the same cardiac-tissue preparation 
can arise because of the interaction of spiral tips with 
small-scale inhomogeneities. 

Conduction inhomogeneities in the ventricle include 
scar tissues, resulting from an infarction, or major blood 
vessels. Some theoretical studies of the effects of tissue 
inhomogeneities have been carried out by using model 
equations for cardiac tissue; however, they have not ad- 
dressed the issues we concentrate on. The interaction 
of an excitation wave with piecewise linear obstacles has 
been studied by Starobin et al 21] to understand the 
role of obstacle curvature in the pinning of such waves. 
Xie et al [2^ 1 have considered spiral waves around a cir- 
cular obstacle and given a plausible connection of the 
ST-RS transition to the size of the obstacle. Panfilov 
et al [13, [13, [1E| have shown that a high concentration 
of randomly distributed non-excitable cells can suppress 
spiral break up. Conduction inhomogeneities can also 
play a very important role in pacing termination of car- 
diac arrhythmias (26[: in particular, it is easier to remove 
a spiral wave once it is pinned to an obstacle, as de- 
scribed in Refs. [13, [1^, than to control a state with 
spiral-turbulence . 

Here we initiate a study that has been designed specif- 
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ically to systematize the effects of conduction inhomo- 
geneities in mathematical models for cardiac arrhyth- 
mias. Our work shows clearly that ST can be suppressed 
or not suppressed by obstacles of different sizes depend- 
ing on where they are placed. As we argue below, this 
sensitive dependence on the sizes and positions of obsta- 
cles must be a manifestation of a complex, fractal-like 
boundary [2^ Is^l between the domains of attraction of 
ST and RS. We also show that inhomogeneities in param- 
eters, which govern ratios of time scales, lead to similar 
results. The models and numerical methods used by us is 
described in section II. Section III contains results; and 
we end with a discussion in Section IV. 



II. MODELS AND NUMERICAL METHODS 

We use the Panfilov ^ [H and Luo-Rudy I [H ^ 
models for cardiac tissue in our studies; the former is 
well suited for extensive numerical studies because of its 
relative simplicity; the latter, being realistic, allows us to 
check that the results we obtain are qualitatively correct 
and not artifacts of the Panfilov model. 

The Panfilov model [sil, [s^l consists of two coupled 
equations, one a partial differential equation (PDE) and 
the other an ordinary differential equation (ODE), that 
specify the spatiotemporal evolution of the scaled trans- 
membrane potential V (denoted by e in Refs. [sil, Is^ ) 
and the recovery variable g, into which this model lumps 
all the effects of the different ion channels: 

dV/dt = VV - f{V) - g; 
dg/dt = eiV,g){kV-g). (1) 

The initiation of action potential is encoded in f(V), 
which is piecewise linear: f{V) = CiV, for V < ei, 
f{V) = -C2V + a, for ei < V < 63, and f{V) = 
C3{V—l),ioTV > 62. The physically appropriate param- 
eters given in Refs. [sT', 32] are ei = 0.0026, 62 = 0.837, 
Ci = 20, C2 = 3, C3 = 15, a = 0.06 and fc = 3. The 
function e{V, g) determines the dynamics of the recovery 
variable: e{V,g) — ei for V < 62, £{V,g) = £2 for V > 62, 
and e{V, g) — £3 for V < ci and g < gi with g-j = 1.8, 
ei = 0.01, €2 = 1.0, and eg = 0.3. As in Refs. 0,^, we 
define dimensioned time T to be 5 ms times dimension- 
less time and 1 spatial unit to be 1 mm. The dimensioned 
value of the conductivity constant £) is 2 cm^/s. 

In spite of its simplicity, relative to the Luo-Rudy 
I (LRI) model described below, the Panfilov model 
has been shown to capture several essential features 
of the spatiotemporal evolution of V in cardiac tissue 
[3i|,[3i,[3l[33|. As in the LR I model the Panfilov model 
also contains an absolute and a relative refractory pe- 
riod. The ways in which spiral patterns appear, propa- 
gate and break up, and the methods by which they can 
be controlled are very similar in these models. To make 
sure that the qualitative features we find are not artifacts 
of the Panfilov model we show explicitly, in illustrative 



cases, that they also occur in the realistic Luo-Rudy I 
model, which is based on the Hodgkin-Huxley formal- 
ism and takes into account the details of 6 ionic currents 
(e.g., Na"*", K"*", and Ca^+) and 9 gate variables for the 
voltage-gated ion channels that regulate the flow of ions 
across the membrane (35| . The concentration difference 
of the ions, inside and outside the cell, induces a poten- 
tial difference of approximately -84 mV across the cell 
membrane in the quiescent state. Stimuli, which raise 
the potential across the cell membrane above -60 mV, 
change the conductivity of the ion channels and yield an 
action potential that lasts typically for about 200 ms. 
Once an action potential is initiated there is a refractory 
period during which the same stimulus cannot lead to 
further excitation. Single cells in the Luo-Rudy model 
are coupled diffusively; thus one must solve a PDE for 
the transmembrane potential V; the time evolution and 
V dependence of the currents in this PDE are given by 7 
coupled ordinary differential equations 35, 36] which we 
give in the Appendix. 

We integrate the Panfilov model PDEs in d spatial 
dimensions by using the forward-Euler method in time 
t, with a time step St — 0.022, and a finite-difference 
method in space, with step size 6x — 0.5 and five-point 
and seven-point stencils, respectively, for the Laplacian 
in cJ=2 and d=3. Our spatial grids consist of square 
or simple-cubic lattices with side L mm, i.e., (2L)'^ grid 
points; we have used L=200. Similarly for the LRI model 
PDEs we use a forward-Euler method for time integra- 
tion, with St = 0.01 ms, a finite-difference method in 
space, with Sx — 0.0225 cm, and a square simulation 
domain with 400 x 400 grid points, i.e., L—9Q mm. We 
have checked in representative simulations on somewhat 
smaller domains that a Crank-Nicholson scheme yields 
results in agreement with the numerical scheme described 
above. 

For both models we use no-flux (Neumann) bound- 
ary conditions on the edges of simulation domain and on 
the boundaries of obstacles. We introduce conduction 
inhomogeneities in the medium by setting the diffusion 
constant D equal to zero in regions with obstacles; in 
all other parts of the simulation domain 13 is a nonzero 
constant. The dimensioned value of D is 2 cm^/s for the 
Panfilov model and between 0.5 cm^/s and 1 cm^/s for 
the LRI model; we use D=0.5 cm^/s in the LRI simula- 
tions we report here . In most of our studies the inhomo- 
geneity is taken to be a square region of side I , with 10 
mm < ^ < 40 mm; however, we have also carried out il- 
lustrative simulations with circular or irregularly shaped 
inhomogeneities. In our three-dimensional simulations 
we use an obstacle of height 4 mm and a square base of 
side 40 mm, i.e., 8 and 80 grid points, respectively (For 
a detailed understanding of the three-dimensional case 
we must also consider the effects of rotational anisotropy 
of muscle fibers in cardiac tissue [l^l, but this lies outside 
the scope of our study.) We also study inhomogeneities in 
which ei in model (1) varies over the simulation domain 
but D is constant. 
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FIG. 1: Panfilov-model spiral turbulence (ST): Transmembrane potentials for two dimension (pseudo-greyscale plots A- 
F) and three dimension (isosurface plots G and H). Two-dimensions:200 mm x200 mm domain and a 40 mm x 40 mm square 
obstacle with left-bottom corner at {x,y). (A) no obstacle -ST; (B) {x = 160, y = 100) ST persists; (c) {x = 150, j/ = 100) 
ST replaced by RS (one rotating anchored spiral); (D) (a; = 140, j/ = 100) spiral moves away (medium quiescent). Three- 
dimensional analogs of (B) and (C): (200 x 200 x 4) domain; an obstacle of height 4 mm and a square base of side 40 mm at 
(E) {x = 140, y = 120, z = 0) and (F) (a; = 140, y = 110, z = 0). 



The initial conditions we use are such that, in the ab- 
sence of inhomogeneities, they lead to a state that dis- 
plays spatiotcmporal chaos and spiral turbulence. For 
the Panfilov model we start with a broken-wavefront ini- 
tial condition: For a system of linear size L at time t=0 
we set g=2, for < x < L and < y < ^, and .9 = 
elsewhere, and y = everywhere except ioi y = ^ + 1 
and < .X < ^, where V ~ 0.9. From this broken wave- 
front a spiral wave develops with a core in the centre of 
the simulation domain and, in the absence of inhomo- 
geneities, evolves to a state with broken spiral waves and 
turbulence (Fig. lA). The spirals continue to break up 
even after 35000 ms for the parameters we use. For the 
LRI model we start from the initial condition shown in 
Fig. 2 A which develops, without an obstacle, into the 
spiral-turbulent state shown in Fig. 2B. 

In the presence of an obstacle the spiral turbulence 

(ST) state of Fig. lA can cither remain in the ST state 
or evolve into a quiescent state (Q) with no spirals or 
the RS state with one rotating spiral anchored at the 
obstacle. Wc explore all these possibilities in the next 
Section. Before we do so, we give the criteria we use to 



decide whether a given state, of the system wc consider, 
is of type ST, RS, or Q. In the Panfilov model, if the 
spiral wave continue to form and break up even up to 
3500 ms, we identify the state as ST (Fig. IB); if, by 
contrast, a single spiral wave anchors to the obstacle and 
rotates around it at least for ten rotations (~ 3500 ms 
for the Panfilov model with a 40 x 40 mm^ obstacle) 
we say that an RS state (Fig. IC) has been achieved (we 
have seen that, once it anchors, this rotation of the spiral 
wave continues even after 100 rotation periods); lastly, if 
the spiral wave moves away from the simulation domain 
and is absorbed at the boundaries within 3500 ms, we 
conclude that the state is Q (Fig. ID). 

For the LRI model, if the spiral formation and break up 
continues upto 2200 ms, we identify the state as ST (Fig. 
2C); if the spiral wave gets anchored to the obstacle and 
completes 4 rotation periods (~ 2200 ms for the obstacle 
we use) we identify the state as RS (Fig. 2D); and we 
say that the state Q (Figs. 2E and 2F) is achieved if 
the spiral wave moves away from the simulation domain 
within 2200 ms. 
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FIG. 2: Luo-Rudy-Model spiral turbulence: Pseudo-greyscale plots in a 90 x 90 mm^ illustrating how the initial condition 
(A) evolves, in the absence of obstacles, to (B) via the generation of spiral waves and their subsequent breakup. In the presence 
of a square obstacle of side I placed with its bottom-left corner at {x,y) we obtain the following: (C)/=18 mm and {x — 58.5 
mm, y — 63 mm) ST persists; (0)1—22.5 mm and {x = 58.5 mm, j/ = 63 mm) RS (one spiral anchored at the obstacle); for 
/=18 and {x — 54 mm,y = 63 mm) spirals disappear leaving the medium quiescent (E) at 800 ms and (F) at 1000 ms. 



III. RESULTS 



Cardiac tissue can have conduction inhomogeneities 
at various length scales. Even minute changes in cell 
or gap-junctional densities might act as conduction 
inhomogeneities(20j : these are of the order of microns. 
Scar tissues or blood vessels can lead to much bigger ob- 
stacles; these are in the mm to cm range so they can be 
studied effectively by using the PDEs mentioned above. 
Here we focus on such large obstacles. As in the ex- 
periments of Ikeda et al we fix the position of the 
obstacle and study spiral- wave dynamics as a function of 
the obstacle size. For this we introduce a square obstacle 
of side I in the two-dimensional (d = 2) Panfilov model 
in a square simulation domain with side L—2QQ mm. We 
find that, with the bottom- left corner of the obstacle at 
the point (50 mm, 100 mm) spiral turbulence (ST) per- 
sists if I < (40 — A) mm, a quiescent state (Q) with no 
spirals is obtained if l~ 40mm, and a state with a single 
rotating spiral (RS) anchored at the obstacle is obtained 
if Z > (40 -|- A) mm. To obtain these results we have 
varied I from 2 to 80 mm in steps of A= 1 mm. Hence 
there is a clear transition from spiral turbulence to stable 



spirals, with these two states separated by a state with 
no spirals. 

The final state of the system depends not just on the 
size of the obstacle but also on how it is placed with 
respect to the tip of the initial wavefront. In our simu- 
lations we find, e.g., that even a small obstacle, placed 
close to the tip [Z=10 mm obstacle placed at (100 mm, 100 
mm)], can prevent the spiral from breaking up, whereas a 
bigger obstacle, placed far away from the tip [1= 75 mm, 
placed at (125 mm, 50 mm)], does not affect the spiral. 

To understand in detail how the position of the obsta- 
cle changes the final state, we now present the results of 
our extensive simulations for the d = 2 Panfilov model in 
a square domain with side 200 mm, i.e., 400 x 400 grid 
points, and with a square obstacle of side l=AO mm. Fig- 
ure 3A shows our simulation domain divided into small 
squares of side Ip mm {lp=10 mm in Fig. 3A). The color 
of each small square indicates the final state of the system 
when the position of the lower-left corner of the obstacle 
coincides with that of the small square: white, black, and 
gray indicate, respectively, ST, RS, and Q. In Figs. 3B 
and 3C we show the rich, fractal-like structure of the in- 
terfaces between the ST, RS, and Q regions by zooming in 
successively on small subdomains encompassing sections 
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FIG. 3: Panfilov-model stability diagram: The effect of an 40 x 40 mm^ obstacle in a 200 x 200 mm'^ domain siiown 
by small squares (side Ip) the colors of which indicate the final state of the system when the position of the bottom- left corner 
of the obstacle coincides with that of the small square (white, black, and gray denote ST, RS, and Q, respectively). (A) for 
lp=10 mm. We get the fractal-like structure of the interfaces between ST, RS, and Q by zooming in on small sub domains 
encompassing parts of these interfaces (white boundaries in A, B, and C with (B) lp=5 mm, (C) lp=2.5 mm, and (D) lp= 1 
mm). 
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FIG. 4: The local time series, interbeat interval IBI, and power spectrum of the transmembrane potential V{x,y,t) at a 
representative point {x,y) in the tissue. When the obstacle is at (160 mm,100 mm) a spiral turbulent state ST is obtained with 
the time series (A), and interbeat interval (B)showing non-periodic chaotic behavior and a broad-band power spectrum (C). 
However, with the bottom-left corner of the obstacle at (150 mm, 100 mm), the spiral wave gets attached to the obstacle after 
9 rotations (~ 1800 ms); this is reflected in the time series (D) and the plot of the interbeat interval(E); after transients the 
latter settles on to a constant value of 363 ms; the power spectrum (F) shows discrete peaks with a fundamental frequency uif 
= 2.74 Hz and its harmonics. Initial transients over the first 50,000 6t were removed before we collected data for calculating 
the power spectrum. 



of these interfaces (white boundaries in Figs. 3 A and 3B) 

and reducing the sizes of the small squares into which we 
divide the subdomain. Clearly very small changes in the 
position of the obstacle can change the state of the system 
from ST to Q or RS, i.e., the spatiotemporal evolution 
of the transmembrane potential depends very sensitively 
on the position of the obstacle. 



The time series of the transmembrane potential 
V{x,y,t) taken from a representative point {x,y) in 
the simulation domain illustrates the changes that oc- 
cur when one moves from the ST to the RS regime in 
Fig. 3. Such time series are shown in Fig. 4. For ex- 
ample, when the obstacle is placed with its bottom-left 
corner at (160 mm, 100 mm), the system is in the spiral- 
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FIG. 5: Inhomogeneities in ei : Inhomogeneities in the parameter ei result in the coexistence of different types of spa- 
tiotemporal behavior in the same system. With e^'^O.Ol and ei"=0.02 (see text), we obtain spatiotemporal chaos outside 
the inhomogeneity but quasiperiodic behavior inside it (A); the latter is illustrated by the power spectrum of V{x,y,t) with 
discrete peaks (B) and the former by a broad-band power spectrum (C). With er*=0.03 and er=0.01 and the left-bottom 
corner of the inhomogeneity placed at {x—140 mm, y—140 mm), single and broken spiral waves coexist in same medium (D), 
whereas, with the inhomogeneity at (a;=60 mm, y—50 mm), a single rotating spiral gets anchored to the inhomogeneity (E, F) 
with quasiperiodic behavior illustrated by the interbeat interval (G) and the power spectrum (H). The power spectrum (H) 
shows six frequencies {uji — 4.06, u!2 = 5.56, LJ3 = 6.57, loa = 7.05, uis = 8.58, and coe ~ 9.07 Hz) not rationally related to each 
other; all other frequencies can be expressed as Yl^^i where the Ui are integers. Initial transients over the first 50,000 3t 

were removed before we collected data for calculating power spectra. 



turbulent state ST. The time series of V from the point 
(51 mm, 50 mm) clearly shows non periodic, chaotic be- 
havior. The times between successive spikes in such time 
series, or interbeat intervals (IBI), are plotted versus the 
integers n, which label the spikes, in Fig. 4B; this also 
shows the chaotic nature of the state ST. Figure 4C shows 
the power spectrum E{uj) of the time series in Fig. 4A; 
the broad-band nature of this power spectrum provides 
additional evidence for the chaotic character of ST. By 
combining Figs. 4A-4C with the pseudo-greyscale plots 
of Figs. lA and IB we conclude that ST is not merely 
chaotic but exhibits spatiotemporal chaos. Indeed, it has 
been shown that the Panfilov model, in the spiral tur- 
bulence regime, has several positive Lyapunov exponents 
whose number increases with the size of the simulation 
domain; consequently the Kaplan- Yorke dimension also 



increases with the system size (see Fig. 4 of Ref.|33j); this 
is a clear indication of spatiotemporal chaos. 

If we change the position of the obstacle slightly and 
move it such that its left-bottom corner is at the position 
(150 mm, 100 mm), the spiral eventually gets attached 
to the obstacle. For this case the analogs of Figs. 4A- 
4C are shown, respectively, in Figs. 4D-4F. From the 
time series of Fig. 4D we see that the transmembrane 
potential displays some transients up to about 2000 ms 
but then it settles into periodic behavior. This is also 
mirrored in the plot of IBI versus n in Fig. 4E in which 
the transients asymptote to a constant value for the IBI 
(363 ms) which is characteristic of periodic spikes. Not 
surprisingly, the corresponding power spectrum in Fig. 
4F consists of discrete spikes at frequencies Um = muj / , 
where m is a positive integer and tj/ is the fundamen- 
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tal frequency (w/ = 2.74 Hz). A simple rotating spiral 
anchored at the obstacle (Fig. IC) will clearly result in 
such a periodic time series in the state RS. 

We do not show the analogs of Figs. 4A-4C for the 
quiescent state Q since the transmembrane potential V 
just goes to zero after an initial period of transients. The 
durations for which the transients last, say in Fig. 4D, 
vary greatly depending on the position of the obstacle 
relative to the spiral tip. We have seen transient times 
ranging from 300 ms to 2000 ms in our simulations. 

We obtain similar results for the three-dimensional 
Panfilov and the two-dimensional Luo-Rudy I models: 
Illustrative pictures from our simulations of spiral turbu- 
lence (ST) and a single rotating spiral (RS) anchored at 
the obstacle are shown in Figs. 1 and 2, respectively. 
From these and similar figures we note that the final 
state, ST, RS, or Q, depends not only on the size of the 
obstacle but also on its position. Obstacles of different 
shapes, e.g., circles, irregular shapes, and two squares 
separated from each other, lead to similar results (see 
www . physics . iisc. ernet . in/~rahul/movies . html for 
representative movies of our simulations). 

We have also explored the effects inhomogeneities in 
parameters such as ei in Eq. (1). In the Panfilov model, 
is the recovery time-constant for large values of g and 
intermediate values of y[32|. As ei increases the absolute 
refractory period of the action potential decreases and 
this in turn decreases the pitch of the spiral wave (cf. 
Fig.3 in Ref.ill]). 

In a homogeneous simulation domain (of size say 200 
X 200 mm^) values of ei > 0.03 produce a single period- 
ically rotating spiral. As ei is lowered, e.g., if ei < 0.02, 
quasiperiodic behavior is seen; this is associated with the 
meandering of the tip of a simple rotating spiral. Even 
lower values of ei, say ei — 0.01 that we have used above, 
lead to spatiotemporal chaos. We now consider an inho- 
mogeneous simulation domain in which all parameters 
in the model except ei remain constant over the whole 
simulation domain. We then introduce a square inhomo- 
geneity inside which ei assumes the value e™ and outside 
which it has the value e'j'"*. Different choices of e™ and 
e™* lead to the interesting behaviors we summarize be- 
low. 

With a square patch of size 40 x 40 mm^, e™ — 0.02, 
and e°"* = 0.01, a spatiotemporally chaotic state is ob- 
tained for most positions of this inhomogeneity. But 
there are certain critical positions of this inhomogene- 
ity for which all spirals are completely eliminated (e.g., 
when the left-bottom corner of the inhomogeneity is at 
a;=70 mm, y=120 mm the spiral moves out of the sim- 
ulation domain). For yet other positions of the inhomo- 
geneity, spatiotemporal chaos is obtained outside the in- 
homogeneity but inside it quasiperiodic behavior is seen 
(Figs. 5A-5C). However, with el" = 0.01 and e?"* = 0.03, 
spiral breakup occurs inside the inhomogeneity and coex- 
ists with unbroken periodic spiral waves outside it (Fig. 
5D), as previously noted by Xie et al [3^. Even in this 
case, for certain positions of the inhomogeneity, a single 



spiral wave gets anchored to it (Figs. 5E, 5F) as in the 
case of a conduction inhomogeneity (Fig. IC). However, 
the temporal evolution of V at a representative point in 
Fig. 5E is richer than it is in Fig. IC: V{x,y,t), with 
x—bl and y=50, displays the interbeat interval of Fig. 
5G; the associated power spectrum shows six fundamen- 
tal frequencies, not rationally related to each other, and 
their combinations; this indicates strong quasiperiodicity 
of V{x, y, t). So, even an inhomogeneity in the excitabil- 
ity of the medium can cause the ST-RS or ST-Q tran- 
sitions we have discussed above for the case of conduc- 
tion inhomogeneities. Furthermore, an inhomogeneity in 
excitability can also lead to rich temporal behaviors as 
shown in Figs. 5 E-H. 



IV. DISCUSSION 

We have shown that spiral turbulence in models of car- 
diac arrhythmias depends sensitively on the size and po- 
sition of inhomogeneities in the medium. In particular, 
we have shown that, with the inhomogeneity at a particu- 
lar position, the state of the spiral wave changes from ST 
to RS as the size of the obstacle increases. We have also 
shown that, for an obstacle with fixed size, this transition 
also depend upon the position of the obstacle. Two im- 
portant questions arise from our work: (1) What causes 
the sensitive dependence of such spiral turbulence on the 
positions and sizes of conduction inhomogeneities? (2) 
What are the implications of our theoretical study for 
cardiac arrhythmias and their control? We discuss both 
these questions below. 

Spiral turbulence (ST) and a single rotating spiral (RS) 
in our models are like VF and VT, respectively, in cardiac 
tissue. Our study suggests, therefore, that such cardiac 
arrhythmias, like their ST and RS analogs in the Pan- 
filov and Luo-Rudy I models, must depend sensitively on 
the positions and sizes of conduction inhomogeneities. 
Furthermore, our work indicates that this is a natural 
consequence of th e sp atiotemporal chaos associated with 
spiral turbulence [33l.[39j in these models: Even for much 
simpler, low-dimensional dynamical systems it is often 
the case that a fractal basin boundary [s^] separates 
the basin of attraction of a strange attractor from the 
basin of attraction of a fixed point or limit cycle; thus 
a small change in the initial condition can lead either to 
chaos, associated with the strange attractor, or to the 
simple dynamical behaviors associated with fixed points 
or limit cycles. 

The PDEs we consider here are infinite-dimensional 
dynamical systems; the complete basin boundaries for 
these arc not easy to determine; however, it is reasonable 
to assume that a complex, fractal-like boundary separate 
the basins of attraction of spatiotemporally chaotic states 
(e.g., ST) and those with simpler behaviors (e.g., RS or 
Q). Here we do not change the initial condition; instead 
we change the dynamical system slightly by moving the 
position, size, or shape of a conduction inhomogeneity. 
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This too affects the long-time evolution of the system as 
sensitively as does a change in the initial conditions. 

In particular, our work elucidates that, by changing the 
position of a conduction inhomogeneity, we may convert 
spiral breakup to single rotating spiral or vice versa as 
depicted graphically in Figs. 3 and 4. Even more exciting 
is the possibility that, at the boundary between these 
two types of behavior (Fig. 3) , we can find the quiescent 
state Q. Thus our model study obtains all the analogs 
of possible qualitative behaviors found in experiments, 
namely, (1) ST might persist even in the presence of an 
obstacle, (2) it might be suppressed partially and become 
RS, or (3) it might be eliminated completely. 

Our work on inhomogeneities in the parameter ei in the 
Panfilov model illustrates the complex way in which the 
spatiotemporal evolution of the transmembrane potential 
depends on the properties of this model for cardiac tissue. 

The implications of our results for anti-tachycardia- 
pacing and defibrillation algorithms, used for the sup- 
pression of cardiac arrhythmias, are very important. Op- 
timal pacing algorithms might well have to be tailor made 
for different inhomogeneities. Indeed, clinicians often 
adapt their hospital procedures for the treatment of ar- 
rhythmias, on a case-by-case basis, to account for cardiac 
structural variations between patients [13]. We hope, 
therefore, that our work will stimulate further system- 
atic experiments on the effects of obstacles on cardiac 
arrhythmias. 
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(plateau K~^), Ii, (total background), given by: 

iNa = GNam^hj{V ~ ENa); 
Is^ - Gs^df{V -E,,); 

Ik = GxxxtiV - Ek); 

Iki — G KiKiooiV — Eki); 

Ikp = GxpKpiV — Ekp); 

h = 0.03921(V^ + 59.87); 

and Kioo is the steady-state value of the gating variable 
Ki. All current densities are in units of fiA/crn^, volt- 
ages are in mV, and G^ and E^ are, respectively, the 
ion-channel conductance and reversal potential for the 
channel ^. The ionic currents are determined by the time- 
dependent ion-channel gating variables h, j, to, d, /, x, 
Xi, Kp and Ki generically denoted by ^, which follow 
ordinary differential equations of the type 



dt 



a^/{a^ + (3^) is the steady-state value of 
is its time constant. The voltage- 



where ^oo 
^ and 

dependent rate constants, and are given by the 
following empirical equations: 



= 0, if y > -40 mV, 
= 0.135 exp[-0.147 (F-l-80)], otherwise; 
1 



if y > -40 mV, 



0.13 (1 -t- exp [-Om{V + 10.66)]) 
= 3.56 exp [0.079 +3.1 x 10^ exp [0.35 V], otherwise; 
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APPENDIX: THE LUO-RUDY MODEL 

In the Luo-Rudy I (LR I) model there are six compo- 
nents of the ionic current, which are formulated mathe- 
matically in terms of Hodgkin-Huxley-type equations [ilj. 
The partial differential equation for the transmembrane 
potential V is 



dV 
'dt 



Ilr 



(A.l) 



Here Ilb. is the instantaneous, total ionic-current den- 
sity. The subscript LR denotes that we use the formu- 
lation of the total ionic current described by the Luo- 
Rudy Phase I (LRl) model 



I si + Ih 



LR -1- J^Kl + Ikp ^ 

(fast inward Na^), Ig 
time- dependent K~^), 



where Ilr = In a + 
/{,, with the current densities I^a 
; (slow inward), Ik (slow outward 
Iki {time-independent K'^)^ Ixp 



0, if y > -40 mV, 

(exp [0.2444 V] + 2.732 x 10-^° exp [-0.04391 V]) 

^ -7.865 X 10-6{1 + exp [0.311 {V + 79.23)]} ^ 

x(T/ + 37.78), otherwise; 

0.3 exp [-2.535 x lO"'^ V] , 

^ ^ \ if F > -40 mV, 



1 + exp [-0.1 (y-f 32)] 
0.1212 exp [-0.01052 V] 
1 + cxp [-0.1378 (y + 40.14)]' 



otherwise; 



ad 



0.32 (T/ + 47.13) 
1- exp [-0.1 (y + 47.13)]' 
0.08 exp [-0.0909 F]; 



0.095 exp [-0.01 {V - 5)] 
1 + exp [-0.072 {V - 5)] ' 
0.07 exp [-0.017 {V + 44)] ^ 
1 + exp [0.05 {V + 44)] ' 



9 



/3/ 



aKi 



Pki = 



0.012 exp [-0.008 {V + 28)] ^ 
1 + exp [0.15 (y + 28)] ' 

0.0065 exp [-0.02 (t/ + 30)]_ 
1 + exp [-0.2 (y + 30)] ' 



0.0005 exp [0.083 {V + 50)] ^ 
1 + exp [0.057 {V + 50)] ' 

0.0013 exp [-0.06 (V + 20)] _ 
1 + exp [-0.04 {V + 20)] ' 

1.02 



1 + exp [0.2385 {V - Eki 
[0.49124 exp [0.08032 {V - 



- 59.215)]' 
Eki + 5.476)] 



1 + exp [-0.5143 {V - Eki + 4.753)] 
-exp [0.06175 {V - Eki - 594.31]]. 



The gating variables Xi and Kp are given by 



2.837 exp0.04(y + 77) - 1 



+ 77) exp 0.04 {V + 
1, otherwise; 

1 

1 + exp [0.1672 (7.488- V)]' 



if F > -lOOmV, 
(A.2) 
(A.3) 



The values of the channel conductances G^a, Gsi, Gk, 
Gk,, and Gkp are 23M).07, 0.705, 0.6047 and 0.0183 
mS/cm^, respectively [42]. The reversal potentials are 
ENa = 54.4 mV, Ek = -77 mV, Eki = Ekp = -87.26 
mV, Eb = -59.87 mV, and Est = 7.7 - 13.0287 In Ca, 
where Ca is the calcium ionic concentration satisfying 



dCa 
dt 



-10"^/sj + 0.07(10"'' -Ca). 



The times t and are in ms; the rate constants ae and 



/3j are in ms 
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